home *** CD-ROM | disk | FTP | other *** search
/ Acorn User: China / Acorn User China CD-ROM (UK) (Disc A) / Acorn User China CD-ROM (UK) (Disc A).bin / HENSA / MATHS / JACAL.ARC / jacal_norm < prev    next >
Encoding:
Text File  |  1993-10-03  |  8.4 KB  |  272 lines
  1. ;;; JACAL: Symbolic Mathematics System.        -*-scheme-*-
  2. ;;; Copyright 1989, 1990, 1991, 1992, 1993 Aubrey Jaffer.
  3. ;;; See the file "COPYING" for terms applying to this program.
  4.  
  5. (proclaim '(optimize (speed 3) (compilation-speed 0)))
  6.  
  7. (define (vsubst new old e)
  8.   (cond ((eq? new old) e)
  9.     ((number? e) e)
  10.     ((bunch? e) (map (lambda (e) (vsubst new old e)) e))
  11.     ((var:> new (car e)) (univ:norm0 new (cdr (promote old e))))
  12.     (else (poly:resultant (make-var-eqn new old) e old))))
  13.  
  14. (define (make-var-eqn new old)
  15.   (if (var:> old new)
  16.       (list old (list new 0 -1) 1)
  17.       (list new (list old 0 -1) 1)))
  18.  
  19. (define (swapvars x y p)
  20.   (vsubst x _$
  21.     (vsubst y x
  22.       (vsubst _$ y p))))
  23.  
  24. ;;; Makes an expression whose value is the variable VAR in the equation
  25. ;;; E or (if E is an expression) E=0
  26. (define (suchthat var e)
  27.   (set! e (poly:subst0 $ (licit->poleqn e)))
  28.   (extize (normalize (vsubst $ var e))))
  29.  
  30. ;; canonicalizers
  31. (define (normalize x)
  32.   (cond ((and math:phases (not (novalue? x)))
  33.      (display-diag 'normalizing:)
  34.      (newline-diag)
  35.      (math:write x *output-grammar*)))
  36.   (let ((ans (normalize1 x)))
  37.     (cond ((and math:phases (not (novalue? x)))
  38.        (display-diag 'yielding:)
  39.        (newline-diag)
  40.        (math:write ans *output-grammar*)))
  41.     ans))
  42. (define (normalize1 x)
  43.   (cond ((bunch? x) (map normalize x))
  44.     ((symbol? x) (eval-error 'normalize-symbol?- x))
  45.     ((eqn? x)
  46.      (poly->eqn (signcan (poly:square-and-num-cont-free
  47.                   (alg:simplify (eqn->poly x))))))
  48.     (else (expr:normalize x))))
  49. (define (expr:normalize p)
  50.   (if (expl? p) (set! p (expl->impl p)))
  51.   (expr:norm-or-signcan
  52.    (poly:square-free-var
  53.     (alg:simplify (if (rat? p) (alg:clear-denoms p) p))
  54.     $)))
  55. (define (extize p)
  56.   (cond ((bunch? p) (eval-error 'cannot-suchthat-a-vector p))
  57.     ((eqn? p) p)
  58.     ((expl? p) p)
  59.     ((rat? p) p)
  60.     (else
  61.      (set! newextstr (sect:next-string newextstr))
  62.      (let ((v (defext (string->var
  63.                (if (clambda? p) (string-append "@" newextstr)
  64.                    newextstr))
  65.             p)))
  66.        (set! var-news (cons v var-news))
  67.        (var->expl v)))))
  68.  
  69. ;(trace normalize normalize1 extize signcan
  70. ;       expr:norm-or-signcan expr:normalize
  71. ;       alg:simplify alg:clear-denoms
  72. ;       poly:square-free-var poly:square-and-num-cont-free)
  73.  
  74. ;; differentials
  75.  
  76. (define (total-diffn p vars)
  77.   (if (null? vars) 0
  78.       (poly:+ (poly:* (var->expl (var:differential (car vars)))
  79.               (poly:diff p (car vars)))
  80.           (total-diffn p (cdr vars)))))
  81.  
  82. (define (chain-rule v vd)
  83.   (if (extrule v)
  84.       (total-chain-exts (total-diffn (extrule v) (poly:vars (extrule v)))
  85.             (var:funcs (extrule v)))
  86.       (let ((functor (sexp->math (car (var:sexp v)))))
  87.     (do ((pos 1 (+ 1 pos))
  88.          (al (cdr (func-arglist v)) (cdr al))
  89.          (sum 0 (app* $1*$2+$3
  90.               (apply deferop
  91.                  (deferop _partial functor pos)
  92.                  (cdr (func-arglist v)))
  93.               (total-differential (car al))
  94.               sum)))
  95.         ((null? al) (vsubst vd $ sum))))))
  96.  
  97. (define (total-chain-exts drule es)
  98.   (if (null? es) drule
  99.       (let ((ed (var:differential (car es))))
  100.     (total-chain-exts
  101.      (poly:resultant drule (chain-rule (car es) ed) ed)
  102.      (if (extrule (car es))
  103.          (union (cdr es) (alg:exts (extrule (car es))))
  104.          (cdr es))))))
  105.  
  106. (define (total-differential a)
  107.   (cond
  108.    ((bunch? a) (map total-differential a))
  109.    ((eqn? a) (poly->eqn
  110.           (total-diffn (eqn->poly a) (poly:vars (eqn->poly a)))))
  111.    (else (let ((aes (chainables a)))
  112.        (if (and (null? aes) (expl? a))
  113.            (total-diffn a (poly:vars a))
  114.            (let ((pa (licit->poleqn a)))
  115.          (total-chain-exts
  116.           (vsubst $ d$ (poly:resultant
  117.                 pa (total-diffn pa (poly:vars pa)) $))
  118.           aes)))))))
  119.  
  120.  
  121. (define (diff a var)
  122.   (cond
  123.    ((bunch? a) (map (lambda (x) (diff x var)) a))
  124.    ((eqn? a) (poly->eqn (diff (eqn->poly a) var)))
  125.    (else (let* ((td (total-differential a))
  126.         (vd (var->expl (var:differential var)))
  127.         (td1 (app* $1/$2 td vd))
  128.         (dpvs '()))
  129.        (poly:for-each-var
  130.         (lambda (v)
  131.           (if (and (not (eq? (car vd) v))
  132.                (var:differential? v))
  133.           (set! dpvs (adjoin v dpvs))))
  134.         td)
  135.        (reduce-init (lambda (e x) (poly:coeff e x 0)) td1 dpvs)))))
  136.  
  137. (define (expls:diff a var)
  138.   (cond
  139.    ((bunch? a) (map (lambda (x) (expl:diff x var)) a))
  140.    ((eqn? a) (poly->eqn (expl:diff (eqn->poly a) var)))
  141.    (else (let ((aes (alg:exts a)))
  142.        (if (and (null? aes) (expl? a)
  143.             (every (lambda (x) (null? (var:depends x))) (poly:vars a)))
  144.            (poly:diff a var)
  145.            (math:error 'polydiff 'not-a-polynomial? a))))))
  146.  
  147. ;(trace total-differential total-chain-exts chain-rule total-diffn
  148. ;       diff expls:diff)
  149.  
  150. ;;;; FINITE DIFFERENCES
  151. ;;; shift needs to go through extensions; which will create new
  152. ;;; extensions (yucc).    It is clear what to do for radicals, but other
  153. ;;; extensions will be hard to link up.  For instance y: {x|x^5+a+b+9+x}
  154. ;;; needs to yield the same number whether a or b is substituted first.
  155. (define (shift p var)
  156.   (vsubst var
  157.       $2
  158.       (poly:resultant (list $2 (list var -1 -1) 1)
  159.               p
  160.               var)))
  161. (define (unsum p var)
  162.   (app* $1-$2 p (shift p (expl->var var))))
  163. (define (poly:fdiffn p vars)
  164.   (if (null? vars) 0
  165.     (poly:+ (poly:* (var->expl (var:finite-differential (car vars)))
  166.             (unsum p (car vars)))
  167.         (poly:fdiffn p (cdr vars)))))
  168. (define (total-finite-differential e)
  169.   (if (bunch? e)
  170.       (map total-finite-differential e)
  171.     (poly:fdiffn e (alg:vars e))))
  172.  
  173. ;;;logical operations on licits
  174. ;(define (impl:not p)
  175. ;  (poly:+ (poly:* (licit->poleqn p)
  176. ;          (var->expl (sexp->var (new-symbol "~")))) -1))
  177.  
  178. ;(define (impl:and p . qs)
  179. ;  (cond ((bunch? p) (impl:and (append p qs)))))
  180.  
  181. (define (expl:t? e) (equal? e expl:t))
  182. (define (ncexpt a pow)
  183.   (cond ((not (or (integer? pow) (expl:t? pow)))
  184.      (deferop _^^ a pow))
  185.     ((eqns? a) (math:error 'Expt-of-equation?:- a))
  186.     ((not (bunch? a)) (fcexpt a pow))
  187.     ((expl:t? pow) (transpose a))
  188.     (else (mtrx:expt a pow))))
  189.  
  190. ;;;; Routines for factoring
  191. (define (poly:diff-coefs el n)
  192.   (if (null? el)
  193.       el
  194.     (cons (poly:* n (car el))
  195.       (poly:diff-coefs (cdr el) (+ 1 n)))))
  196. (define (poly:diff p var)
  197.   (cond ((number? p) 0)
  198.     ((eq? (car p) var) (univ:norm0 var (poly:diff-coefs (cddr p) 1)))
  199.     ((var:> var (car p)) 0)
  200.     (else (univ:norm0 (car p) (map-no-end-0s
  201.                    (lambda (x) (poly:diff x var))
  202.                    (cdr p))))))
  203. (define (poly:diff-all p)
  204.   (let ((ans 0))
  205.     (do ((vars (poly:vars p) (cdr vars)))
  206.     ((null? vars) ans)
  207.     (set! ans (poly:+ (poly:diff p (car vars)) ans)))))
  208.  
  209. ;;;functions involved with square free.
  210. (define (poly:square-free-var p var)
  211.   (poly:/ p (poly:gcd p (poly:diff p var))))
  212.  
  213. (define (poly:square-and-num-cont-free p)
  214.   (if (number? p) (if (zero? p) p 1)
  215.       (poly:* (poly:square-and-num-cont-free (univ:cont p))
  216.           (poly:/ p (poly:gcd p (poly:diff p (car p)))))))
  217.  
  218. (define (poly:factorq p) (poly:factor-all p))
  219.  
  220. (define (poly:factor-var c var)
  221.   (poly:factor-split c (poly:diff c var)))
  222.  
  223. (define (poly:factor-all c)
  224.   (poly:factor-split c (poly:diff-all c)))
  225.  
  226. ;;; This algorithm is due to:
  227. ;;; Multivariate Polynomial Factorization
  228. ;;; David R. Musser
  229. ;;; Journal of the Association for Computing Machinery
  230. ;;; Vol 22, No. 2, April 1975
  231.  
  232. (define (poly:factor-split c splitter)
  233.   (let ((d '()) (aj '()) (b (poly:gcd c splitter))) ; changed #f's to ()
  234.     (do ((b b (poly:/ b d))
  235.          (a (poly:/ c b) d))
  236.         ((number? b)
  237.          (if (eqv? 1 b)
  238.              (cons a aj)                ; nreverse removed
  239.              (cons b (cons a aj))))     ; nreverse removed
  240.       (set! d (poly:gcd a b))
  241.       (set! a (poly:/ a d))             ; 'a' is used as a temporary variable
  242.       (if (not (eqv? a 1))              ; keeps extraneous `1's
  243.         (set! aj (cons a aj))))))       ;  out of the results list
  244.  
  245. ;;; the following algorithm attempts to separate factors in a multivariate
  246. ;;; polynomial with major variable.  It substitues 0 for each variable
  247. ;;; that it finds in turn and takes GCD against the original expression.
  248. ;;; It assumes that it's argument is squarefree and contentfree in the
  249. ;;; major variable.
  250. (define (univ:split pe varlist)
  251.   (cond ((unit? pe) (list))
  252.     ((null? varlist) (list pe))
  253.     ((let ((p0 (signcan
  254.             (poly:gcd pe (poly:subst0 (car varlist) pe))))
  255.            (cvl (cdr varlist)))
  256.        (if (unit? p0)
  257.            (univ:split pe cvl)
  258.          (nconc (univ:split (poly:/ pe p0) cvl)
  259.             (univ:split p0 cvl)))))))
  260.  
  261. (define (univ:split-all poly) (univ:split poly (poly:vars poly)))
  262.  
  263. (define (factor-free-var p var)
  264.   (poly:gcd p (poly:subst0 var p)))
  265.  
  266. (define (factor:test)
  267.   (test (list (list x -1 -2) (list x -1 0 1))
  268.         poly:factorq
  269.         (list x -1 -4 -3 4 4)))
  270.  
  271. ;;;    Copyright 1989, 1990, 1991, 1992, 1993 Aubrey Jaffer.
  272.